Nonlinear r-Modes in Neutron Stars: Instability of an unstable mode 
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We study the dynamical evolution of a large amplitude r-mode by numerical simulations. i?-modes 
in neutron stars are unstable growing modes, driven by gravitational radiation reaction. In these 
simulations, r-modes of amplitude unity or above are destroyed by a catastrophic decay: A large 
amplitude r-mode gradually leaks energy into other fluid modes, which in turn act nonlinearly with 
the r-mode, leading to the onset of the rapid decay. As a result the r-mode suddenly breaks down 
into a differentially rotating configuration. The catastrophic decay does not appear to be related 
to shock waves at the star's surface. The limit it imposes on the r-mode amplitude is significantly 
smaller than that suggested by previous fully nonlinear numerical simulations. 



PACS numbers: 95.30.Sf, 04.30.Db, 04.40.Dg, 97.60.Jd 



The r-modes of rotating neutron stars are unstable 
growing modes driven by gravitational radiation reaction 
pi 0- If t ne ^ — rn = 2 r-mode of a young, rapidly ro- 
tating star can grow to an amplitude of order unity, the 
gravitational radiation it emits would carry away most of 
the star's angular momentum and rotational kinetic en- 
ergy; and the radiation might be detectable by the Laser 
Interferometric Gravitational Wave Observatory (LIGO 
II) Q (for recent reviews see 0,0). Even if its am- 
plitude were smaller, the r-mode instability would limit 
the periods of hot, young neutron stars (and, possibly, 
of old stars spun up by accretion). A number of mech- 
anisms to damp the mode have been examined, includ- 
ing shear viscosity enhanced by crust-core coupling and 
by nonstandard cooling 6]; bulk viscosity enhanced by 
a hyperon-rich core ; and energy loss to a magnetic 
field driven by differential rotation 8] . But none of these 
definitively eliminates the instability. 

The significance of the r-mode instability depends 
strongly on its maximum possible amplitude. In two 
recent Letters, Stergioulas and Font [!j and Lindblom, 
Tohline, and Vallisneri |10| performed numerical simu- 
lations of nonlinear r-modes, both finding that large- 
amplitude nonlinear r-modes can exist for some long pe- 
riod of time. In addition, in |lOt ITU ]. Lindblom et al. 
carried out numerical simulations of the growth of the r- 
modes driven by the current quadrupole post-Newtonian 
radiation reaction force in Newtonian hydrodynamics. In 
order to achieve a significant growth of the r-mode am- 



plitude in a reasonably short computational time, they 
artificially multiplied the radiation reaction force by a 
factor of 4500. This decreases the growth time of the 
r-mode from about 40 s to 10 ms. The (dimensionless) 
r-mode amplitude a grew to ~ 3.3 before shock waves 
appeared on the surface of the star and the r-mode am- 
plitude collapsed. Lindblom et al. suggested that the 
nonlinear saturation amplitude of the r-modes may be 
set by dissipation of energy in the production of shock 
waves. Here, we show that a hydrodynamical effect will 
restrict the r-mode amplitude to a value significantly be- 
low that reported in [nj, [llj . 

We note that, with the artificially large radiation reac- 
tion, the results in |l0llll| assume that no hydrodynamic 
process takes energy from the r-mode in a time scale be- 
tween 10 ms and 40 s (the artificial growth time and the 
actual physical growth time, respectively). In this paper 
we investigate the evolution of large amplitude r-modes 
in these time scales (10 ms-40 s). We find that (i) a 
catastrophic decay of the mode, occuring within these 
time scales, significantly reduces the amplitude to which 
an r-mode can grow, and (ii) in a large amplitude r-mode 
this catastrophic decay leads to a differentially rotating 
configuration. (As in Refs. jTnL Hll | , we assume a perfect 
fluid with no magnetic fields.) 

We solve the Newtonian hydrodynamics equations for 
a non-viscous fluid flow in the presence of gravitational 
radiation reaction: 
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P (J + g ' W ) = ~ VP ~ PV$ + pi?GR ' (2) 
where p is the density, P is the pressure, v is the velocity, 
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Fgr is the radiation reaction force per unit mass, and 
is the Newtonian potential, satisfying 



V 2 $ = AirGp 



(3) 



A high resolution shock capturing scheme (Roe solver) 
is used to solve the hydrodynamic equations. In addition, 
as in |9j, we ap plie d the 3rd order piecewise parabolic 
method (PPM) 12] for the cell-reconstruction process in 
order to simulate rapidly rotating stars accurately for a 
large number of rotational periods. 

As r-modes couple to the gravitational radiation 
mainly through the current multipoles J; m , we assume 
that the contribution to the reaction force -Pgr comes 
solely from the dominant current multipole J22, an ap- 
proximation used also in |lOt . The resulting expres- 
sion for F GR is given by (see [13111111) 



F, 



GR 



GR 



fr GR — 



-Kl (x + iy) 



3v z J. 



(5) 



z J, 



(6)' 



k Im < (x + iyY 



-J. 



(5) 



r(6) 



>y 



(n) 

where J 22 is the nth time derivative of J22 and k — kq = 
32y / 7rG/(45v / 5c 7 ). A technical difficulty in evaluating 
-Fgr is that it depends on high-order time derivatives of 

(n) 

J 22- To circumvent this problem, we assume that J 22 = 
(iw)" J22, with the nonlinear r-mode frequency oj defined 
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be expressed as an integral over the fluid variables and is 
thus calculated on each time slice [ljOj, [llj, ll4j . 

To investigate the dynamical properties of a large am- 
plitude r-mode numerically, one must first generate ini- 
tial data. However, analytically we know only the ex- 
pression of an r-mode in the linear regime, under the 
assumption of a <C 1 and £! — > 0. In particular, the 
/ = 77i — 2 r-mode perturbation (at the first order of fig) 
is given by 



(6) 



where «o is a dimensionless amplitude, Rg and f2 are the 
radius and angular velocity of the unperturbed rotating 
star model, and Y 2 % is the magnetic-type vector spher- 
ical harmonic defined by = [1(1 + l)] _1 'Vx VY; m . 
Using Eq. © with a large ciq introduces other modes in 
the initial data |Tlj : that is, the resulting configuration 
is not the same as an r-mode growing from small ampli- 
tude driven by gravitational radiation reaction. In the 
following, we will show how we obtain a large amplitude 
r-mode in our study. 

Beyond the linear regime, there is no unique mode 
decomposition, and hence no unique definition of an 
r-mode. In this paper, by a large amplitude r-mode 
we mean a configuration resulting from the growth of 




FIG. 1: Evolution of a in a slowly rotating star with the 
correct (ko) and artificial (9 x 10 7 ko) radiation reaction. The 
solid lines represent the numerical results (257 3 resolution), 
while the dashed lines are the predictions from linear theory. 



an infinitesimal r-mode due to gravitational radiation 
reaction. We define the amplitude a of the nonlin- 
ear I = m = 2 mode in terms of its contribution to 
J22 = / pr 2 7j-Y 2 B 2 *d 3 x: 



*(*) 



a(x)e l4 ' {x) = 



8nR* (pr 2 v-Y£*) 
n(t) J pr 4 d 3 x 



(7) 



Here R is the radius of the corresponding nonrotating 
star model and a(x) is the amplitude density. The phase 
factor 4>(x) is defined so that a(x) is real. The definition 
of the amplitude is similar to that of except that 
we normalize a(t) by the average angular velocity of the 
star O(i) instead of a fixed initial f2o- 

In Fig. ^ we show the time evolution of a for two 
cases in a slowly rotating star with T = 4.42 ms. In the 
first case k is set to the correct Post-Newtonian value kq 
(represented by the solid line labeled "ko"). The evo- 
lution begins with a small (linear) r-mode perturbation 
given by Eq. © with ao = 0.1. The simulation is car- 
ried out up to t — 0.6 ms. To compare, we plotted as 
dashed line the evolution of a as predicted by the lin- 
ear theory a — aoe t/TGR , where the gravitational 
radiation time scale is given by 
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We note that this formula is correct only to linear order in 
a and in Q, and is thus accurate only for (i) a << 1 and 
(ii) flo/^/Gp << 1. While assumption (i) is reasonably 
good, (ii) is actually not accurate for the model used: 
rotation period 4.42 ms, corresponding to f2/f2 max « 0.25 
(^max ~ 2^/ttGp/S is the Kepler limit, where p represents 
the average density of the star). Nevertheless, the two 
lines nearly coincide. 

With radiation reaction coefficient kq, 
time t ps 3 x 10 s ms is needed to reach 1 



129 grid-point simulation, this would take O(10 ) time 



an evolution 
t = 1. For a 
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FIG. 2: Left: Growth of the r-mode amplitude a during arti- 
ficial "pumping" with k = 4500ko- Right: The velocity profile 
v y along the x axis at t = and at the point when a = 2.0. 



steps, requiring a clearly impractical O(10 8 ) hours on a 
128 CPU Origin 2000 (MIPS R12000). 

To arrive at a large am plit ude r-mode we use an ar- 
tificially large k as in [1(1 Hlj. In Fig. ^ the solid line 
labeled "9 x 10 7 Ko" represents the evolution of a with 
K = 9 x 10 7 Ko- For comparison, we also plot the evolu- 
tion of a as predicted by the linear theory with this k 
as the dashed line. The slight offset between the dashed 
and solid lines at the initial time is due to the fact that 
a(t = 0) as calculated from Eq. (J7J) equals to ao as de- 
fined in Eq. © only in the limit a ^ 1 and ~~ * 0. 

In Fig. [21 (left), we show the growth of an r-mode 
(with ao = 0.5) to a large amplitude with an artificial 
k of 4500ko m our fast rotating star model: The star 
has a mass of M = 1.64M Q with a poly tropic equation 
of state P = kp 2 . The equatorial radius R e — 14.5 km. 
The ratio of the polar to equatorial radii is 0.76. The 
rotation period is T = 1.24 ms. This model is used for 
the rest of the simulations discussed in this paper. Unless 
otherwise noted, we use 129 3 Cartesian grid points with 
Ax = Ay = Az = 0.42 km. 

In Fig. [21 (left) , the amplitude a rises from 0.5 to 2.2 in 

19 ms. An indicator of the accuracy of the simulation is 
that the total mass of the system is constant to 0.08% by 

20 ms in our 129 3 runs. Also, the actual numerical evo- 
lution of the total angular momentum J = \J pfx vd 3 x\ 
agrees with the theoretical prediction [see Eq. (11) of 
Ref. [13] to about 1%. 

While the discussion above shows the accuracy of our 
numerical treatment, one still must ask whether the 
large amplitude r-mode obtained with the large artificial 
pumping is physical or not, in the sense that whether the 
rapid pumping excites modes that would not be excited 
with k — kq . We compare the large amplitude r- modes 
obtained with different pumping rates and conclude that 
the resulting fluid flow pattern does not depend sensi- 
tively on the pump rate (as long as the pump rate is large 
enough so that the large amplitude r-mode can be arrived 
at). In Fig. [21 (right), we show the rotational velocity pro- 
file v v along the x axis for n = 9000, 4500ko at the point 
when a = 2.0, starting with the same initial model. For 
comparison, we also plot the initial profile in the same fig- 
ure. We see that the two lines, k — 9000, 4500«o, agree 
with each other to better than 3%, with smaller discrep- 



FIG. 3: Left: Evolution of the "pumped up" a at different 
values with k — kq. Right: Fourier spectra of v z along the x 
axis (at x — 6 km) at two time slots in the evolution starting 
out with a = 1.6. The early (later) time slot is denoted by 
the dashed (solid) line. 



ancy away from the surface. In the rest of this paper, we 
will take the "pumped up" configurations given in Fig. [21 
(left) as the initial state in our investigation of the hydro- 
dynamical behavior of the large amplitude r-mode. To 
the extent that different pumping rates are not affect- 
ing the initial state we use, the artificial pumping is not 
affecting the results we report below. 

Here we focus on one question: What is the fate of a 
large-amplitude r-mode in a rapidly rotating neutron star 
for a sufficiently long time evolution, modeled as the 1.24 
ms polytrope described above. That is, what evolution 
is implied by Eqs. 0-1®, with the correct amount of 
radiation reaction? 

In Fig. [21 (left) we show the evolution of a vs time for 
various large amplitude r-modes starting off with a = 
2.2, 2.0, 1.8, 1.6. We see that the mode amplitudes 
start off slowly decaying, leaking energy to other modes. 
The decay rate is small, until a certain time. We plot 
in Fig. [21 (right) the Fourier transform of the velocity 
component v z along the x axis at a typical point inside 
the star at two different time slots in the evolution of the 
a = 1.6 case [the lowest line in Fig. [21 (left)]. The dashed 
line is the profile at the beginning of the evolution (13 
ms-35 ms), we see that there is only one large peak at 
the r-mode frequency (0.93 kHz). This is compared to 
the solid line representing the spectrum at a later time 
slot (35 ms - 50 ms). We see that various smaller peaks 
appear in the spectrum, especially the one at twice the 
r-mode frequency (1.86 kHz). Spectra at different points 
inside the star give similar structure, with those in the 
core region showing more peaks at different frequencies. 

The most interesting feature of Fig. [21 (left) is that 
after some slow leaking of energy into other fluid modes, 
the r-mode amplitude drops catastrophically to a value 
much smaller than 1. This abrupt drop occurs through 
nonlinear couplings with other fluid modes: In the slow 
leaking phase, these other fluid modes are growing lin- 
early until a certain unstable point. The time it takes 
to reach the unstable point depends sensitively on the r- 
mode amplitude. It shortens from approximately 45 ms 
to 8 ms when the initial value of a changes from 1.6 to 
2.2. 
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FIG. 4: Evolution of a (65 3 resolution) with artificial pump- 
ing of k = 4500kq whenever a drops below its initial value. 



To further investigate this catastrophic decay we per- 
form a set of numerical experiments in which we pump 
energy into the r-mode by turning on the artificial radi- 
ation reaction force with coefficient k = 4500ko when- 
ever its amplitude drops down below its initial value. 
The resulting evolution of a vs time is given in Fig. 
0] The evolution tracks of a starting off with values 
a = 2.1, 1.9, 1.7, 1.5, 1.4, 1.0 are given. With the 
large artificial pumping a remains constant despite en- 
ergy leaking to other modes. In all cases (except a = 1.0 
where the simulation is not evolved long enough), how- 
ever, the hydrodynamical nonlinear interaction eventu- 
ally overwhelms the artificial pumping, and the r-mode 
amplitude falls catastrophically. 

In Fig. [S]we compare the distribution of the amplitude 
density a(x) as defined in Eq. JJJ) on the equatorial plane 
before (left) and after (right) the breakdown respectively 
for the case where the initial a is 2.0 (the a = 2.0 line 
in Fig. [3J|. In the figure, the brighter region represents 
higher amplitude density. During the catastrophic de- 
cay, the r-mode pattern changes rapidly from a 4-fold 
"regular" shape (left) to a whirlpool- like spiral (right). 
We also see in our simulations that strong differential ro- 
tation is developed during the breakdown, a potentially 
important fact regarding whether subsequent re-growth 
of the r-mode is possible or not |16( . 

In Fig. E|we plot the rotational velocity profile v v along 




FIG. 5: The amplitude density a(x) on the equatorial plane 
before (left) and after (right) the breakdown for the case 
where the initial "pumped up" a is 2.0 in Fig. |H] 



FIG. 6: The rotational velocity profile v y along the x axis for 
the case where the initial "pumped up" a is 2.0. The solid 
line is the initial profile at a = 2.0, while the dashed line is 
the profile after the breakdown. 



the x axis for the case where the initial a is 2.0. The solid 
line is the initial profile, while the dashed line is the pro- 
file after the decay. To further quantify the amount of dif- 
ferential rotation, we define the kinetic energy associated 
to differential rotation by / = \ J p (v^ — v^) cPx, where 

= Cly^x 2 + y 2 with Cl being the average angular veloc- 
ity of the star. We plot the evolution of a and / together 
in Fig. [7| (left). It is seen that the amount of differential 
rotation (I) rises rapidly during the breakdown of a. We 
also see in our simulations that the star has a relatively 
large amplitude pulsation during the breakdown. Fig. 
{7\ (right) shows the quadrupolc-moment component Q xy 
against time for the same case as in Fig. (left). We 
see that Q xy is basically zero until the breakdown and 
it then oscillates rapidly afterward. However, based on 
an order-of-magnitude estimation, the gravitational radi- 
ation amplitude due to the changing quadrupole moment 
is only about 1% of that due to the r-mode. 

In contrast to the study of Refs. 0, [n|, we do not 
see evidence that this catastrophic decay is due to the 
generation of shock waves on the surface of the star. 

We found in this paper that a large amplitude r-mode 
will lose energy to other fluid modes whose growth in 
turn trigger a catastrophic decay of the r-mode. For an 
r-mode with an amplitude a of order one, the onset of 
catastrophic decay requires a time much shorter than the 
growth time of the r-mode due to radiation reaction; and 
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FIG. 7: Left: Evolution of a (dashed) and I (solid) for the 
case where the initial "pumped up" a is 2.0. Note that I has 
been rescaled for comparison. Right: Evolution of Q xy (in 
G — c — Mq = 1 units) for the same case. 
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the decay thus limits the r-mode amplitude to a value less 
than that found bv [ToLIT]|. Further work is in progress to 
determine the nature of this catastrophic decay, whether 
it is related to any known hydrodynamical instability of 
nonlinear flow, and how large an r-mode amplitude can 
be with this effect taken into account. 

Note added. Towards the end of the preparation of 
this paper we learned the results of Arras et al. 01 • We 
note the following differences between the hydrodynami- 
cal phenomenon studied in this paper and that studied by 
them: The work of Arras et al. points to a slow leakage 
of the r-mode energy into some short wavelength oscil- 
lation modes, leading to an equilibrium distribution of 
mode amplitudes. This in turn, through viscosity dissi- 
pation, limits the r-mode amplitude to a small value. In 
this paper we find a sudden and complete breakdown of 
the r-mode that operates independent of viscosity. Fur- 
ther numerical investigation will be carried out to inves- 
tigate the interactions of r-modes with other oscillation 
modes. Such investigation is beyond the resolution power 



of our present simulations. Note that, if the conclusions 
of are correct, there may be no astrophysical situa- 
tion in which r-modes grow to amplitudes large enough 
to exhibit the sudden decay seen in our simulations. 
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